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Abstract 

In this contribution a numerical study of a turbulent jet flow is presented. The simula- 
tion results of two different variants of the Lattice Boltzmann method (LBM) are com- 
pared. The first is the well-established D3Q19 MRT model extended by a Smagorinsky 
Large Eddy Simulation (LES) model. The second is the D3Q21 Factorized Cascaded 
Lattice Boltzmann (FCLB) model without any additional explicit turbulence model. 
For this model no studies of turbulent flow with high resolution on nonuniform grids 
existed so far. The underlying computational procedure uses a time nested refinement 
technique and a grid with more than a billion DOF. The simulations were conducted 
with the parallel multi physics solver VirtualFluids. It is shown that both mod- 
els are feasible for the present flow case, but the FCLB outperforms the traditional 
approach in some aspects. 

Keywords: lattice Boltzmann, large eddy simulation, factorized cascaded Lattice 
Boltzmann, distributed simulation, jet 



1. Introduction 

The Lattice Boltzmann model can be considered an alternative approach to obtain 
numerical solutions of the Navier-Stokes equations, even though LBM can also be used 
to investigate finite Knudsen number flows. LBM is based directly on the distribution 
functions for the particle dynamics of the fluid. The method has successfully been 
employed to model and simulate a variety of complex fluid flow problems ranging 
from multi component [1] and multi phase flows [2] to thermal flows [3], fluid-structure 
interaction [4], non-Newtonian flows [5] and turbulent flows [6]. Over the last years 
a a number of Lattice Boltzmann variants have been developed to simulate turbulent 
flows [6,7]. 

Even though Direct Numerical Simulation (DNS) is gaining more relevance for 
certain turbulence flow problems, it is still prohibitively expensive for most relevant 
applications including turbulence. Any mature CFD scheme should also be capable 
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of incorporating state-of-the-art turbulence models. In the Lattice Boltzmann context 
large eddy simulation (LES) models are particularly popular due to the small timestep 
of the explicit scheme and the small overhead needed to implement an algebraic LES 
model [7, 8], but RANS models have also been used with LBM [6]. 

An alternative approach to the simulation of turbulent flows using turbulence mod- 
els is the use of numerical methods without any explicit turbulence model but relying 
entirely on a suitable dissipation of the numerical scheme. The fine turbulent scales 
are not resolved, and the numerical discretization is acting as a filter Such schemes, 
named implicit large eddy simulation (ILLS) models, are becoming more popular as 
stated by Grinstein et al. [9]. The Factorized Cascaded Lattice Boltzmann (FCLB) 
model has been shown to give reasonable results at high Reynolds numbers with very 
low resolution [10] without any explicit turbulence model. 

Our simulations are based on the research code VirtualFluids - a parallel code 
which is based on MPI and the METIS partitioning tool [11]. A hybrid block data 
structure to overcome the bottlenecks of the previous approach is used [12]. This block 
data structure enables partitioning of very large datasets because only the block data 
structure has to be partitioned instead of the entire set of individual nodes. For local 
grid refinement with hierarchical block grids [12] the grid refinement strategy of Yu et 
al. [13] is employed. See also [14] for a review and evaluation of various refinement 
techniques and Crouse et al. [15] for applications. 

Jet flow is a standard validation problem that has been studied thoroughly both 
experimentally and numerically, such as in the early experimental work of Wygnanski 
et al. [16] and DNS study of Boersma et al. [17]. A Lattice Boltzmann study of a 
turbulent square jet flow has been carried out by Yu et al.[18, 19]. The MRT (eq. (12)) 
and SRT (eq (1)) models with Smagorinsky LES have been compared on a uniform 
grid with a D3Q19 stencil, a 19-element stencil in three directions. [20] conducted a 
further study of a square jet with Lattice Boltzmann LES. 

In this article two different Lattice Boltzmann collision models, namely the D3Q19 
MRT model with Smagorinsky LES and the D3Qn Factorized Cascaded Lattice Boltz- 
mann (FCLB) model, are evaluated for their capability to predict turbulent flows for the 
complex flow case of a free jet. For axisymetrical flows a lack of isotropy has been re- 
ported for the D3Q15 and D3Q19 models, while the D3Q21 was found to remove this 
flaw as White and Chong [21] observed when they tested the isotropy of these lattices 
for flow through a nozzle at Re < 500 using the BGK and MRT model. They pointed 
out the importance of reducing isotropy errors as they had found that the errors de- 
pended only weakly on the grid resolution. Mayer and Hazi [22] also observed a lack 
of isotropy for the D3Q19 but not the D3Q21 model in a study of laminar and turbulent 
flow through rod bundles. 

The article is structured as follows: We start with an overview over different LBM 
variants, the D3Q19 MRT model, the D3Qn models Cascaded Lattice Boltzmann 
(CLB) and FCLB. The incorporation of large eddy models in LBM is briefly recalled. 
In the second part of the article we present the testcase of the turbulent jet flow. Firstly, 
the flow type, for which a semi-analytical solution is known, is described. Next we 
give a brief description of the experiment to which we compare our data. After that the 
numerical setup is presented followed by simulation results for the D3Q\9 MRT model 
with Smagorinsky LES and the D3Q11 FCLB model. Finally, the results are discussed 
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and differences between the results from the two approaches are pointed out. 



2. Lattice Boltzmann collision models and subgrid stress model 

The Lattice Boltzmann scheme emerged in the late 1980's from Lattice Gas Cellu- 
lar Automata [23] as a new approach to Computational Fluid Mechanics. Unlike con- 
ventional discretizations of the Navier Stokes equations. Lattice Boltzmann equations 
rely on a discretization of a simplified Boltzmann equation which is a time-dependent 
description of the behavior of particle ensembles. In its simplest form, it is based on a 
single relaxation time for the non-equilibrium distribution function [24]. 

Mx + eiAt,t+At)-fi=^{Mx,t)-f^\x,t)) (1) 

for a distribution function /, its equilibrium f'' and the relaxation parameter T. The 
components fi{x,t) of the distribution function depend on the discrete time step t, the 
position X which is related to a discrete node of the numerical grid and the index ; for 
the discrete velocity set. The equilibrium for the incompressible model reads [25] 

f,'> = ^,^5p+3-^ + -—, -^j (2) 

Here 5p is the density fluctuation for p = po + 5p, u is the macroscopic velocity and 
Cs the speed of sound in the LBM context. By e we denote the discretized microscopic 
velocity. For the D3Q19 model the weight factors are wq = 1 /3, vv'i = 1/18, vv2 1 /36 
and for the D3Q27 model wq = 8/27, wi = 2/27, W2 = 1/54, W3 = 1/216 where W3 
is used for the velocity vectors that point to the corners of the cube. The entries of the 
velocity vectors e, for the D3Q19 and D3Q27 model are: 
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where c is the lattice speed At /Ax. The macroscopic values density p and velocity 
u are computed from summation of the distributions. For the incompressible model we 
set p = po + 5p and 

Sp = L/, (3) 

Ux = —y\fieix (4) 
Po i 

Uy = -^Jlf'^iy (5) 
POi 

for incompressible models (such as the D3Q19 MRT LES model we used in this work) 
or 

P = Lf- (6) 

Pux = J^Mix (7) 

PUy = 52 //e;,. (8) 

for compressible models (such as the D3Q27 FCLB model used in this work). More 
advanced approaches use Multiple Relaxation Times (MRT) [26] where the relaxation 
step takes place in moment space. To introduce moments let us first define an expec- 
tation value of a linear operator B acting on distribution functions f in the discretized 
velocity space 

w = L(fi(/)),- (9) 

Moments are then defined as expectation values of powers of the discrete velocities 

IX^^j^, - (e'^e^et^ (10) 

in accordance with the definitions for distribution functions in continuous spaces that 
can be found in e.g. [27]. An alternative notation that avoids multiple subscripts is 

Mj^^2^3...3-=^.vVz* (11) 

'■ J k 

The density and momentum defined in eqns. (6) to (8) are the moments of order zero 
and one. The transformation from distribution functions to moments is a linear trans- 
formation M. The Lattice Boltzmann relaxation step for the MRT model is then given 
by the equation 

fi{x + ei,t + At)= fi+M-'^S{Mf{x,t)^nf''{x,t)) (12) 

with nf^ = Mfi. The collision parameters s„ of the diagonal matrix S are chosen 
such that the density and momentum are conserved and that the viscosity is correctly 
represented using 

V 1 

T = 3^ + -Ar, (13) 
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where V is the kinematic viscosity in lattice units [v] = j ISt. Additional free pa- 
rameters can be chosen as to improve the stability of the model. The moments, or 
central moments, belong to different invariant subsets of the stencils' symmetry groups 
as described in reference [28]. They have to be relaxed with one relaxation factor each 
(definition see below). Table 1 lists these groups and the corresponding relaxation fac- 
tors. For the D3219 model only the moments for s\, S2, .?3, S4, and 57 are considered. 
We chose the values of the relaxation parameters to be i,- = 1 V/ 7^ 1 and si = Af /t for 
both the D3Q19 and the D3Q21 simulation runs. 



Table 1 : relaxation factors 



moment 


relaxation parameter 


^12, Mn, M23, Ml 1 - M22, Ml 1 - M33 


si 


Mil +M22+M33 
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Ml23 
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MII22 — 2Mi133 +M2233, Ml 122 + Ml 133 ^ 2M2233 


Si 


MII22 + M1I33 +M2233 


Si 


M1223. M1223, M1233 


Si 


M12233. Ml 1233. Ml 1223 


S9 


Ml 12233 


sw 



The first order moments and the density do not appear in table 1 as they are con- 
served. 

A further development are the cascaded Lattice Boltzmann schemes. The so-called 
Cascaded Lattice Boltzmann (CLE) method was developed by Geier et al. [29]. Fur- 
ther developments have been made to introduce different equilibria [10] which consti- 
tute the FCLB scheme. All CLB-methods rely on the basic idea to use central moments 
instead of uncentered moments and to use lower order moments after relaxation for the 
computation of the higher order moments (hence the term cascaded). Central moments 
are defined as 

M^; = ((x-(x))') (14) 

for the expectation value (x) of a function f{x). In our case, the expectation value 
is intended for the discrete distribution function f{x,Jl,t) with respect to momentum 
space as defined above. For three directions we have a product of one-dimensional 
terms. 

= {{e,-{e,)Y{ey^{e,)yie,-{e,))') (15) 

Intuition leads us to suspect that non-linear operations should be carried out in the iner- 
tial frame. Consider for example the second order central moments (i.e. the variances). 
The uncentered moment is Mi 1 = + var{vT(). Hence the term Y^ift^'^ix^ix ^ 9 1^ 
has been removed by the transformation and the central moment then is the variance 
only. The equilibrium central moments are chosen as the corresponding central mo- 
ments of the Gauss function where the variance is the speed of sound Cj. These are 
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the same equilibria as those obtained from taking the central moments of the MRT- 
equilibria if third-order terms are taken into account for the MRT equilibria as well. 
The Factorized CLB method is a special CLB method which aims at removing the in- 
fluence of the lower-order central moments on the fourth- and higher order moments 
at an acceptable computational cost. This correction leads to an improved stability of 
the method and further reduces errors with respect to isotropy that occur with any finite 
stencil [10]. The transformation and specific equilibria for the D3Qn stencil are given 
in table 2. The original implementation of Geier et al. [10, 29] computed the changes 
in the moments after collision. Our implementation differs from the original imple- 
mentation as we do not compute the change in the moments, but recompute the entire 
moments. The basis for the moments used in reference [29] has some differences from 
the basis used here. 

We chose this implementation because of its more modular properties. The first 
transformation is the same as for the MRT model. The less compressed implementa- 
tion is less prone to errors and makes it easier to change algorithmic details later. On 
the other hand, it is not as optimized as the original version with respect to the number 
of floating point operations (FLOPS). A large number of FLOPS can be eliminated 
if relaxation parameters are fixed. The CLB and FCLB model are suspected to have 
ILES capabilities. This has been subject to investigation in references [29], [10] and 
[30] where no additional turbulence model was used. For under-resolved simulations of 
turbulent flows with the LBGK or MRT model, however, a turbulence model is needed. 
The standard Smagorinsky model is a popular choice due to its simplicity and effi- 
ciency. In this model the eddy viscosity Vx depends only on the magnitude of the strain 
rate S and the grid spacing Ax 

Vt = (CsAx)2||S|| (16) 



where the strain rate tensor is defined as 



1 / dua du 



lydxp^ dxa J ' ^^'^^ 

and the Smagorinsky constant Cs- We chose Cs = 0.18, which is in the range of values 
suggested by Rogallo and Moin [31]. In the Lattice Boltzmann context the viscosity v 
is related to the relaxation time t as defined in equation (13). The norm of the strain 
rate can be computed locally from 

llsiH-^r-^lin-^ll, (18) 

^ HotalC 

where the norm of the momentum flux tensor n"''i is defined as 

K', -(e(C;-5p/35„,)') ' (19) 

in the case of an incompressible model. The total relaxation factor can be obtained 
from the following equation [8] 



Ttotal = ^Vq + -M+ ^ ^ (20) 
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Table 2: transformation to central moments 
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where 

Note that the procedure is entkely local. No information from adjacent nodes is re- 
quired, which is highly desirable for parallel computations. For the description of the 
hierarchical block structured grid approach for the SGS model we refer to reference 
[32]. 



3. Validation of turbulent jet flow 

A turbulent jet at Re = 6760 based on the size of the orifice and the inflow velocity 
is simulated using the FCLB method with the D3Q27 stencil and with the MRT model 
with Smagorinsky LES and the D3Q19 stencil. The simulation results are compared to 
experimental data from Ming et al. [33]. The section is structured as follows: Firstly, 
the experimental setup is described. The setup of the numerical solution is described 
after that, followed by the results of the simulations. Finally, the results are discussed 
and differences between the results from the two approaches are pointed out. 

3.1. Experimental setup 

The simulations are based on an experiment described in reference [33]. The prop- 
erties of a turbulent jet at a Reynolds number of 6760 based on the size of the opening 
of 4mm and on the inflow velocity of 1.69m/5 were measured using Doppler laser 
anemometry. The experiment was carried out in a water tank of 6m length in flow 
direction, 0.2m width and 0.4m height. The tank is open and the jet enters the tank 
through a nozzle. At the back of the water tank a drain is present to keep the water 
level constant. 



3.2. Numerical setup 

With the numerical setup we try to mimic the experimental setup as closely as pos- 
sible. We use the same size of domain in horizontal, vertical, and spanwise direction. 
Solid boundaries are modeled by no-slip boundaries. The air-water interface at the 
upper boundary is modeled by a free-slip condition because a free-surface condition 
would pose a major additional computational effort and the effect of the wave genera- 
tion is considered to be negligible for this testcase. Instead of the weir outflow we set a 
fixed pressure boundary condition. The nozzle was positioned at 0.5m, approximated 
as a cylinder with second-order accurate interpolated no-slip walls [34]. The point of 
origin is on the bottom, left, frontal corner of the basin. Instead of the nozzle used in the 
experiment a cylinder is inserted, which extends from (0.0, 0.1, 0.2) to (0.5, 0.1, 0.2) 
meters and has a radius of 2mm. On the right emitting end of the cylinder a constant 
inflow velocity is defined. The boundary condition at the walls is a noslip condition. 

For the discretization of the domain a hybrid block structured grid with a hierarchi- 
cal refinement structure is used. Due to the geometrical refinement a nested time step 
approach is used leading to a globally constant CFL number for the distributions. The 
refinement and coarsening strategy is described in [4, 12, 35]. Seven levels of refine- 
ment are used to discretize this setup. The grid resolution is 0.0947mm on the finest and 
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6.06051mm on the coarsest level. So the nozzle with its 4mm of diameter is discretized 
with 42.24 nodes in the finest domain. The timestep varies between 0.000103522i 
(coarse) and 0.0000016^ (fine). The domain is resolved with 83522 blocks, each of 
which corresponds to nodal matrix of the size 11x11x11. In sum 111 million grid 
nodes are used. This means three billion degrees of freedom for the D3Q21 model and 
2. 1 billion for the D3Q19 model. The domain was decomposed for parallelization with 
the METIS library [36]. 
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Figure 1: partial view of the domain discretization witli bloclcs of 11x1 1x1 1 nodes each , the color indicates 
the subdomain index after decomposition with METIS 

The physical parameters are the fluid density of 998.2 kg/nv', the kinematic viscos- 
ity of 10^^ m^/.s, the Reynolds number (Re) 6760 (related to the nozzle diameter and 
inflow speed), and the computation time which covered 3.9i real time. 

3.3. Results 

We compare the averaged velocity along the axis abtained for the FCLB and for 
the MRT model with the semi-analytical results from [33]. Figure 5 shows a good 
match for both models. Pictures 4 and 7 give a qualitative idea of the flow dynamics. 
Immediately behind the opening the flow field is laminar As eddies develop in the 
shear layer between jet and surrounding flow, the jet becomes wider with increasing 
distance from the nozzle. 

According to Ming et al. [33] the average axial velocity behind the nozzle can be 
described as: 



D 



(22) 



with nozzle diameter D = 4mm, the distance x from the nozzle, the averaged ve- 
locity Um at position x, and the inflow speed uq — l.69m/s. The constant k\\ has to be 
determined experimentally and was determined to kn = 6.104 for the present setup. 



Figure 2: system: pipe with velocity contour at 0.5m/s after time t = 3.9s 




average vx 1 
1..9439a6 




-0.31112 



Figure 4: averaged horizontal velocity behind the nozzle, D3Q19 MRT LES 




Figure 5: averaged velocity along the jet axis 
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The spreading width bg for the jet is defined as the half-width of the velocity over 
the distance from the jet axis at a given distance from the jet assuming a Gaussian 
shape. This leads to a velocity = 0.368m,„ which is present at half the spread-width 
from the jet axis. The spreading width grows linearly with the distance from the jet 
[33] 

bg = kj^x (23) 

Ming et al. [33] found a value of /tx = 0.109. 

Due to the asymmetric behavior of the jet, the minimum and maximum radius of the 
spreading function for the distance to the isoline of constant velocity u^^ at = 0.368m„, 
is given in figure 6. 





Figure 7: averaged velocity ortliogonal to jet axis 5cm behind nozzle, left D3Q19 MRT LES, right D3Q27 
FCLB, contour at 0.368Mm 

The velocity distribution orthogonal to the jet axis is determined according to [33] 

by: 

«, = M„,- 0.938 (24) 

with the averaged velocity at jet axis m,„, the velocity Ux at the position r, the ra- 
dius re where = 0.368m,„, and radius r. The constants have again been determined 
experimentally. For different distances behind the nozzle figures 9, 10 and 1 1 show the 
computed results in comparison with the semi-analytical solution. For the computation 



12 




Figure 8: averaged velocity profile 1.5cm behind nozzle, left D3Q19 MRT LES, right D3Q27 FCLB 



of the averaged velocities as well as the turbulent intensity, several lines in different 
directions from the jet center orthogonal to the jet axis are averaged in addition to 
averaging in time. In figure 12 the distribution of the turbulent intensity at different 
positions is shown which was determined from 





Figure 9: averaged velocity profile 5cm behind nozzle, left D3Q19 MRT LES, right D3Q27 FCLB 




Figure 10: averaged velocity profile 9cm behind nozzle, left D3Q19 MRT LES, right D3Q27 FCLB 



Figures 5 to 1 1 show the computed results for the two models in comparison with 
the semi-analytical solution. As can be seen from figure 5, the D3Q27 FCLB model is 
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Figure 11: averaged velocity profile 14a« behind nozzle, left D3Q19 MRT LES, right D3Q27 FCLB 




Figure 12: turbulent intensity orthogonal to jet axis, left D3Q19 MRT LES. right D3Q27 FCLB 
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Figure 13: experimental data for the turbulent intensity, from [33] 
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slightly more successful at reproducing the velocity profile along the jet centerline than 
the D3Q19 MRT model with Smagorinsky LES. The same is true for the spreading 
width (fig. 6) for moderately large distances from the nozzle. For distances larger 
than 7cm the error in the spreading width of the FCLB model grows, but this may 
be due to the limited averaging time of 0.9 seconds. The same behavior is observed 
for the average velocity profiles normal to the jet axis. We believe that the excessive 
eddy viscosity that occurs with the Smagorinsky LES model in shear layers delays the 
transition to turbulence. 

An interesting observation is that the mean velocity contours normal to the jet axis 
diverge from the expected circular shape for the D3Q19 model. The discretization of 
the velocity space with 19 vectors seems insufficient to reproduce this particular flow 
feature. The use of the D3Q27 FCLB model improves the isotropy of the flow field. 
Similar effects have been observed previously by White and Chong [21] in a compar- 
ison of D3Q19 and D3Q21 BGK-type models at Reynolds numbers up to Re — 500. 
[10] shows a comparison between different stencils and collision models for laminar 
flows and also found that the D3Q21 FCLB model showed the least anisotropy among 
the models studied. 

4. Conclusion 

In this paper we presented a comparison of a D3Q19 MRT model with Smagorin- 
sky LES and the D3Q21 FCLB model. We demonstrated that both models correctly 
reproduce the dynamics of turbulent jet flow. The computation of one second real time 
on 395 cores took two days. The decay of the axial velocity is in good agreement 
with the semi-analytical solution. The solution from D3Q27 FCLB model matches the 
semi-analytical result better than the D3Q19 LES model. In the range of 0cm to 10cm 
behind the nozzle the spreading functions are in good agreement with the empirical 
relation determined from experiments. The velocity profile of a cross-section matches 
the Gauss function obtain from empirical relations well. One important aspect is that 
the D3Q19 LES model shows notable anisotropics whereas the D3Q27 FCLB model 
shows no such defect. 

We conclude that the Lattice Boltzmann method is suitable for jet induced turbulent 
incompressible flows even with a simple turbulence model (LES) and an enhanced 
model (FCLB) used in this work. The potential of the FCLB model for computing 
turbulent flows is demonstrated. 
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